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Abstract. We investigate the spreading of information in a one-dimensional Bose- 
Hubbard system after a sudden parameter change. In particular, we study the time- 
evolution of correlations and entanglement following a quench. The investigated 
quantities show a light-cone like evolution, i.e. the spreading with a finite velocity 
We discuss the relation of this veloctiy to other characteristic velocities of the system, 
like the sound velocity. The entanglement is investigated using two different measures, 
the von-Neuman entropy and mutual information. Whereas the von-Neumann entropy 
grows rapidly with time the mutual information between two small sub-systems can 
as well decrease after an initial increase. Additionally we show that the static von 
Ncuman entropy characterises the location of the quantum phase transition. 
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1. Introduction 

Entanglement and correlations are important properties of quantum systems. The 
interest in these quantities is manifold. Correlations have been used since a long time to 
characterize the low energy properties of quantum many body states. In particular in the 
study of quantum phase transitions, correlations played a leading role. Entanglement 
has been the resource for quantum computing and quantum cryptography. Only recently 
entanglement has also become a major tool to understand and characterize ground states 
of interacting quantum many body systems. Driven by these applications the study of 
entanglement has become an active field of research and many different measures have 
been proposed in quantum many body systems [1]. 

Recently the interest in the fundamental questions of quantum dynamics was 
reinforced by the experimental realization of strongly interacting gases [2]. In these 
systems the precise and rapid tunability of the system parameters and the very good 
decoupling from the environment open the possibility to study the quantum evolution of 
a system far from equilibrium. In this context the following questions arise: How does 
a quantum system react to a parameter change? How fast can information propagate 
through the system? How do correlations and entanglement between different parts of 
the system build up? 

These and similar questions have been asked in different situations. Lieb and 
Robinson found for one-dimensional spin systems that the quantum dynamics generated 
by local operators obeys a light-cone like evolution, i.e. the information spreads with a 
finite velocity [3]. They proved that only exponentially small corrections exist outside 
this light-cone in the considered systems. Many generalizations of this Lieb Robinson 
theorem have been developed over the last years (see Ref. [4] and references therein). 
Further the finding of the spreading of correlations with a finite velocity after a global 
quench has been verified in specific integrable models [5, 6, 7, 8, 9], as well as for critical 
systems which can be described by conformal field theory [10, 11]. 

Motivated by the experimental investigation of a sudden parameter change across 
the superfluid to Mott-insulator transition [12] we investigate here the time-evolution 
of a one-dimensional non-integrable Bose-Hubbard system after a sudden quench of the 
interaction strength. For this model the recent generalization of the Lieb-Robinson 
theorem of Cramer et al. [4] applies deep in the superfluid regime. However, no 
generalization exists in the presence of sizable interaction strength between the bosons. 

We start our investigation with the static properties of the von Neumann 
entanglement entropy in the Bose-Hubbard model. In particular we show that previous 
predictions from conformal field theory [13, 14] are in very good agreement with our 
numerical results in the critical superfluid phase assuming a central charge c — 1. We 
further point out how the knowledge on the dependence of the entanglement entropy 
on the system and block length in the critical phase can be used to locate the quantum 
phase transition. 

We then turn to the discussion of the time-evolution of the quantum system after 
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a sudden parameter change. Hereby we mainly consider a quench from the superfluid 
to the Mott-insulating phase. However, we also show some examples of quenches inside 
the superfluid or the Mott-insulating regime. In a previous work [15] we focussed on the 
relaxation of the system to a quasi-stationary state after a quench. Here in contrast we 
discuss the short time behaviour in detail. We study the spreading of information by the 
time-evolution of correlations, the von Neumann entropy, and the mutual information. 
We find that for onsite quantities the time-scale of relaxation is set by the hopping time 
of the bosons. In contrast for longer range correlations, like the one-particle density 
matrix or the density-density correlations, a front spreads with almost constant velocity 
which causes a 'light-cone' like evolution. The front travels at a finite speed, such that 
in an infinite system equilibrium can never be reached. We discuss the time evolution of 
the von Neumann entropy of different contiguous blocks of length I. We find that the von 
Neumann entropy of a certain block shows a sharp linear rise for short time and saturates 
for longer times. To a good approximation the rate of the block entropy increase is 
shown to follow a boundary law, while the entropy value at saturation depends only on 
the block volume. Our results nicely corroborate recent analytical results [16, 17, 18]. 
In contrast to the von Neumann entanglement, the mutual information between two 
contiguous blocks of different length after a quench to a Mott-insulating regime drops 
substantially in time compared to the values in the initial superfluid state. So for the 
considered quench the growth of the entropy and the loss of correlations go hand in 
hand. 

2. Model and methods 

Ultracold bosons subjected to an optical lattice can be described by the Bose-Hubbard 
model [19, 20, 21], here in its one-dimensional form 

H(J, U) = -JJ2 (b]b j+1 + h.c) + ! 5>,(n, - 1), (1) 

j j 

where &t and bj are the boson creation and annihilation operators, and rij = bjbj 
the number operators on site j. The first term in Equation (1) models the kinetic 
energy of the atoms and the parameter J is called the hopping parameter. The second 
term stems from the short range interaction between the atoms and the parameter U 
characterizes its strength. In the experimental setup of ultracold bosons in an optical 
lattice the parameters U and J can be directly related to the experimental parameters. 
In particular, by varying the lattice height in the experiment the parameter U / J can be 
changed by several orders of magnitude and by anisotropic lattices the dimensionality of 
the system can be varied. In equilibrium at integer filling the Bose-Hubbard model shows 
a quantum phase transition at a critical value of u c := {U / J) c between a superfluid state 
(U/J < u c ), in which the atoms are delocalized, to a Mott insulating state (U/J > u c ) 
in which the atoms are localized [21]. We use h — 1 and a — 1, where a denotes 
the lattice spacing, throughout the work. The theoretical tools we use to treat the 
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ground state and the dynamic properties of the Bose-Hubbard model are analytical 
approximations and numerical methods. Static results for the entanglement entropy are 
obtained using the density matrix renormalization group method (DMRG) [22, 23, 24] 
with up to 1000 states. The results for time-evolution are calculated using the exact 
diagonalization (ED) based on Krylov methods [25, 26] and the adaptive time-dependent 
DMRG (t-DMRG) [27, 28, 29]. The ED is used to study one-dimensional chains with 
periodic boundary conditions up to 16 sites, while the adaptive t-DMRG is used for 
one dimensional systems with open boundaries and with up to 64 sites keeping several 
hundred DMRG states. 

3. Static von Neumann Entanglement Entropy 

2.2 
2 

1.8 
1.6 

1.4 | 
1.2 c/) -1 
1 

0.8 
0.6 

500 1000 1 2 3 4 5 6 7 

| Log Conformal Distance 

Figure 1. Left panel: The static von Neumann block entropy Sl(1) for an open 
system of length L = 1024 and a range of different interaction values U/J, located 
in the superfiuid and the Mott insulating phases. The average filling is n = 1. 
Right panel: the same data as a function of the logarithmic conformal distance 
A := log [(2L/ir) sin(7rZ/L)]. The linear behavior of the curves of U/J = 1, 2 and 
3 reveals the c = 1 critical theory. The rapid saturation of the entropy for U/J = A 
and 5 is a consequence of the short correlation length in the Mott insulating phase. 

Before we discuss the behaviour of the system after a sudden change of the 
parameters we would like to discuss the static properties of entanglement in the Bose- 
Hubbard model. In recent years many different measures have been proposed for the 
entanglement in a system [1]. One of the proposed measures is the von Neumann block 
entropy. Let A be a block of length I and B be the remaining system. Then the von 
Neumann entropy of a block A is defined by Sa = — Tr^ pA log pa, where pA = Tr# p is 
the reduced density matrix of the block A and p = \^>){^>\ the pure-state density matrix 
of the whole system. Tr^ denotes the trace over X. 

In the superfiuid phase the low energy physics of the one-dimensional Bose-Hubbard 
model can be described by a Luttinger liquid, which is a conformal field theory with 
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central charge c = 1. For such a 1+1 critical system the von Neumann entropy has 
been derived for different geometries of the subsystems A and B [13, 14]. If A is a single 
block of length I in a periodic system of length L the von Neumann entropy is given by 



Here s\ is a non-universal constant. In a system with open boundary conditions which 
is divided at some interior point into an interval A of length / and its complement the 
von Neumann entropy is 



Here log g is the boundary entropy of Affleck and Ludwig [30] . Note that an oscillating 
correction term beyond the conformal field theory prediction (3) has been found for 
open boundary conditions in the critical phase of the 5 = 1/2 XXZ spin model [31]. 
This is a specific example of Friedel-like oscillations decaying away from the boundaries. 
These, in turn, lead to an oscillating, algebraically decaying correction term to the block 
entropy. In the case of the Bose Hubbard model considered here, which is not particle- 
hole symmetric, the open boundaries can induce a slightly non-uniform particle density, 
but the amplitude of the difference with respect to the nominal density decays rapidly 
away from the boundaries. The effect is most pronounced for small U / J and becomes 
negligible for large U/J. This inhomogeneous particle density distribution is therefore 
expected to affect the entanglement entropy on small open systems. 

For stronger interactions the quantum phase transition to the Mott-insulating 
phase takes places. The Mott-insulating state is characterized by an energy gap above 
the ground state and has a finite correlation length £. For gapped systems with a 
finite correlation length £, the block entropy is expected to saturate at a finite value 



We show in Fig. 1 the von Neumann entropy of the static Bose-Hubbard model 
with open boundary conditions for different interaction strengths U/J obtained using 
DMRG calculations on a system of length L = 1024. In the left panel we show the 
entropy as a function of the block length I, while in the right panel we rescale the a;- axis 
according to the logarithmic conformal distance A := log [(2L/n) sin(7rZ/L)]. Eq. 3 then 
simplifies as Sz,(A) = c/6 A + logg + si/2. 

For small interaction strength U/J = 1,2,3 our results agree very well with the 
prediction (3) of the conformal field theory for large block length I. For U/J = 1 
deviations at small I can be seen. They are due to the slightly inhomogeneous density 
in the open system. Linear fits to the curves in the right panel for the larger values of A 
yield slopes which are very close to the expected value of 1/6. The extracted value of c for 
the three critical U / J values is shown close to the corresponding curve in the left panel. 
The values of c are very accurately 1. The constant contribution logg + s\/2 becomes 
smaller for increasing U/ J due to the reduced onsite number fluctuations. For U / J = 
4 and 5 a very clear saturation of the entropy as a function of the block length I can 




(2) 




(3) 



S L (l) ~ log(£) for I » £ [32, 14]. 
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be seen, revealing the finite correlation length in the Mott insulator. The entanglement 
entropy vanishes completely in the limit of U/ J — ► oo when the groundstate becomes a 
trivial product wave function of singly occupied sites. 
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Figure 2. Difference AS(L) = Sl(L/2) - S L / 2 {L/4) versus U/J for different system 
lengths L at average filling n — 1. For U/J < u c the expected scaling behaviour 
AS = c/61og2 with c = 1 is seen. In the L — »• oo limit the data converges to a 
step function located at u c . The vertical line shows the previously obtained critical 
values u c ~ 3.37(12) [33] (black line with grey uncertainty) and u c m 3.380(4) [34] (red 
line) for comparison. The deviation of small system length for U < 2J stems from 
inhomogcncous density distribution due to the open boundaries and is not captured in 
the CFT approach. Inset: AS as a function of 1/L for selected values of U/J close to 
the phase transition. 



In the following we investigate whether it is possible to locate the critical value u c 
purely based on the properties of the von Neumann entropy. Earlier high precision work 
with DMRG mostly used a priori knowledge about the Kosterlitz-Thouless transition 
to accurately locate the critical value [33, 34]. A recent investigation based on quantum 
information related quantities such as the fidelity and single site entropies and their 
derivatives arrived at the conclusion that the single site entropy would not allow to 
locate the critical value u c precisely [35] . 

In order to locate the critical point we devise the quantity 

AS(L) := S L (L/2) - S L/2 {L/A) (4) 

the increase of the entropy at the mid-system interface upon doubling the system size. 
In a region of parameter space described by a conformal field theory with central charge 
c one expects simply AS* = c/6 log 2 based on Equation (3). On the other hand in a 
gapped region with a finite correlation length one obtains AS* = for L ^> £, because 
the entropy saturates for block lengths I larger than the correlation length £. As the 
system size increases one therefore expects AS for the one-dimensional Bose-Hubbard 
model to scale to a step function as a function of U / J. 
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In Figure 2 we display AS as a function of U/J for different system sizes L from 
32 up to 1024 sites. In the Luttinger liquid regime at small U/J the curves for different 
system sizes nicely converge towards the expected value of 1/6 log 2. For larger U/J 
values the AS" curves scale to zero for increasing L, indicative of the gapped phase. 
Based on our available data we can safely infer that the point U/J = 3.4 is already in 
the Mott insulating phase, while U/J = 3.3 is still critical, see inset of Figure 2. This 
result is in full agreement with the previously obtained critical values u c ~ 3.37(12) [33] 
and u c ~ 3.380(4) [34], based on fits of the decaying bosonic Green's function. In future 
studies using a fine grid of U/J values one could perform a finite size scaling of the 
inflection point to extract an even more accurate value of u c . So we conclude that an 
appropriate scaling plot of the von Neumann entanglement entropy provides competitive 
results on the location of the quantum phase transition in the one-dimensional Bose 
Hubbard model without relying on a priori knowledge of the Kosterlitz-Thouless nature 
of the transition. 

4. Description of the parameter quench 

We implement the quench by an abrupt change of the interaction strength from an 
initial value U to a final value Uf at fixed hopping parameter J at time t = 0. In most 
cases we start from a superfluid phase (U/J < u c ) and change to the Mott-insulating 
regime (Uf/J > u c ). The initial wave function \ip Q ) for t < is the ground state of the 
Hamiltonian Hi = H(J, U/). We investigate its time evolution for times t > subject to 
the Hamiltonian Hf = H(J,Uf). The time evolution of the wave function is governed 
by the Schrodinger equation, i.e. 



The time evolution of expectation values of relevant operators can be expressed as 



Hereby \m) are the eigenstates of the final Hamiltonian Hf and E m the corresponding 
energy eigenvalues. The expansion of the initial state |V>o) into the eigenstates of the 
final Hamiltonian, i.e. \ip ) = ^2 m c m \m) with the weights c m , is used. In a realistic 
case the tunneling between lattice sites and the interaction strength are non-zero and to 
determine the time-evolution of the wave function is not an easy task. We calculated the 
time-evolution numerically using ED and DMRG methods and analytical approximation 
in certain limits. 

5. Light cone effect in correlations 

In this section we investigate how correlations over a distance r react to the sudden 
parameter change. We consider two different types of correlations, the single particle 



#(*)) = exp(-i# / *)|*o)- 




(5) 




Figure 3. Time-evolution of correlation functions after a quench from Ut = 2J to 
Uf = 40J. The upper panel shows the single particle correlation functions (&o& r ) 
for different distances r. The correlations show partial revivals up to a time t r when 
they start to reach a quasi-steady state. This time t r grows approximately linearly 
with the distance r as marked by the vertical lines. The central panel shows the 
same correlations functions after filtering out the high frequencies, see text for details. 
The lowest panel shows the density density correlations function (non r } after shifting 
and rcscaling their amplitude for better visibility. The common vertical dashed lines 
denote the arrival of the minima as determined from the density-density correlations. 
The data shown is ED for a L — 14 and DMRG data for L = 32 and filling n = l. 



correlations (bjbj +r ) and the density- density correlations (rijUj+r) at equal time. In 
Fig. 3 we show the time-evolution of the different correlations after a quench from the 
superfluid, Ui = 2, to the Mott-insulating, Uf = 40, parameter regime. 

Single-particle correlations The upper panel shows the correlations {b\b r ) for different 
distances r\. For short times the single particle correlations oscillate with the period 
27r/Uf. The origin of these oscillations lies in the integer spectrum of the operator 
hj(hj — l)/2. Consider the limit of very strong interactions, where the time-evolution 
is totally dominated by the interactions. The time evolution of the single particle 
correlations is given by 

{m},{m'} 

Here we use the notation {m} for the Fock state with raj particles on site i. The time- 
evolution of the correlation function is determined by the non-vanishing cross terms 

| To extract these correlations from the DMRG data with open boundary conditions the average over 
central sites is taken. Note that for periodic boundary conditions this quantitiy is real due to symmetry, 
whereas for open boundary conditions an imaginary part can develop. However for the shown functions 
and times the imaginary part is negligible. 



Correlations and entanglement after a quench in the Bose-Hubbard model 



9 




Frequency co [J] Frequency co [J] 

Figure 4. Fourier transform of the equal-time single-particle and density-density 
correlation functions for different distances r. The same parameters as in Figure 3 
are used. In the single-particle correlation function clear frequency bands located 
at multiples of the interaction strength Uf — 40 J can be seen. The density-density 
correlation is dominated by the contributions at low frequencies from zero up to ~ J. 

({m}\blbj\{m'}) of Fock states whose occupations vary by removing one particle from 
site j and adding it at site %. The frequencies occuring in the time-evolution will be 
given by U fivn!- — vnl i — 1). This results in oscillations of period T « 2k /Uf for equally 
occupied lattice sites j and i of the state \{m'}). For unequal onsite occupations higher 
multiples of Uf can occur in the frequency spectrum. 

Thus the distribution of the different frequencies in the Fourier transformation of the 
single particle correlation functions (left panel in Fig. 4) and the occupation difference 
in the initial state are intimately connected. For the shown quench ([/, = 2 J, Uf = 40 J) 
sizable contributions of the lowest three frequency bands can be seen, whereas the 
occupation of higher frequency bands becomes very small. The distribution among the 
frequency bands changes by varying the initial value of the interaction strength U%. If 
the initial state is deep in the superfluid regime, the peaks at higher frequencies show 
more weight due to the presence of strong particle number fluctuations. In contrast if 
the initial state is close to the Mott-insulating transition the particle fluctuations are 
suppressed and the weight of the peaks at higher frequencies decreases accordingly. The 
width of the frequency bands is due to the finite value of J and is responsible for the 
decay of the oscillations in time [15]. 

Let us now come back to the real-time evolution of the single particle correlations 
for different distances. After a time t « 0.5/ J the correlation corresponding to the 
smallest distance shown, r = 2, reaches a quasi-steady value (marked by the leftmost 
vertical line). The correlations of distance r = 3 deviate from the correlation with r > 3 
for a time t ~ 0.7/ J (second left vertical line). The same can be observed for correlations 
with increasing distances at longer times. In the central part of the Figure 3 the same 
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single particle correlations are shown, but the high frequency oscillations at frequencies 
u n ~ n x Uf with n> 1 are filtered out, and only the n = contributions are kept. In 
these low-pass filtered correlations the saturation to a quasi-steady state value can be 
seen much more clearly. 
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Figure 5. Time-evolution of the rescaled density-density correlations (non r )(t) after 
a quench with the same parameters as in Figure 3. The front of the evolution evolves 
approximately with a constant velocity. 



Density- density correlations The lowest panel of Figure 3 shows the density- density 
correlations (njUj +r ) — (nj)(nj +r ) at equal time. The amplitudes of the correlations are 
rescaled and shifted for better readability §. 

The density-density correlations do not show strong oscillations, but remain almost 
constant in time up to the moment, where a pronounced signal arrives. In their Fourier 
spectrum (right panel of Figure 4) mostly low frequencies ~ J occur. This is due to 
the fact that the interaction term of the Hamiltonian commutes with the correlation 
function. Thus, in the strong coupling limit the interaction term does not give rise to 
oscillations. Approximately at the same time as the single particle correlations saturate, 
the spreading of a signal (here the reaching of a minimum, loci of the common vertical 
dashed lines in Figure 3) can be found in the density-density correlation (non r ). 

In Figure 5 we show a contour plot of the rescaled density correlations. In this 
representation a clear light cone evolution can be seen, i.e. a front travels through the 
system at almost constant speed. Let us note, that the same light cone effect occurs in 

§ We will denote the rescaled function (njjij +r ^ — (%)(^j+r) averaged over different sites j in the 
center of the chain by (ngn,.). 
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the evolution of the correlations for different quench parameters (cf. Table 1) and also 
in incommensurate systems, e.g. at filling n = 1/2. 
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Figure 6. Characteristic time i r for arrival of the signal in the density-density 
correlations depending on the distance r are summarized for different quench parameter 
Ui, Uf, and average density n. A clear linear relation between time and distance is 
observed for all shown quenches. To give an order of the uncertainties, typical error 
bars are plotted for chosen points. They take into account both the difficulty to identify 
a sharp signal and deviations between different system length (L — 14 up to L — 64). 



Propagation velocity v s To extract the velocity v s of the signal propagation in the 
density-density correlations, we plot in Figure 6 the time t r , when a signal, e.g. the 
minimum, occurs in the density- density correlations versus the distance r. The curves 
for different values of the initial and final interaction strength Ui and Uf are presented. 
A clear linear behaviour of the time t r versus the distance is seen. The shift between 
the curves, e.g. for Uf = 20,40 and Uf = 4,6, stems from the different signatures that 
have been tracked in the density-density correlations. The inverse of the slope is the 
signal propagation velocity v s of the signature. We extract the velocity v s by a linear 
fit t r = r/v s + b, where v s and b are the fitting parameters. The results for v s are shown 
in Table 1. 

The velocity of the spreading of correlations after a quench has been identified 
by Calabrese and Cardy [10] within conformal field theory and for different integrable 
models to be twice the maximal mode velocity. The simple picture given is that the 
modes depart from both considered sites and the signal arrives, if both modes interfere. 
In the description by a conformal field theory this velocity agrees with the sound velocity 
of the system, since the dispersion relation is linear. In a chain of harmonic oscillators 
or in a spin chain lattice effects occur, which cause a curved dispersion relation. Thus 
the maximal velocity of a mode can be distinct from the sound velocity in the system. 
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Figure 7. Dependence of the velocity v s /2 on the final interaction strength after 
the quench. Results are shown for different values of the density and compared to an 
analytical approximation described in the text. The initial value is Ui = 2 J. 



Table 1. Velocity extracted from a linear fit. 
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For the transverse Ising model a velocity equal to one has been found [5] . 

In Figure 7 we show the results for the rescaled signal velocity v s /2 for one initial 
interaction strength U{ = 2.1 . We present results for different densities. At the 
commensurate density n = 1 a phase transition to the Mott-insulator takes place at 
a critical value. In contrast for the incommensurate density n = 1/2 || the system stays 
in the superfluid phase for all interaction strengths. The signal velocity shows a strong 
increase for low interaction strength. After reaching a maximum around U/J m 6 it 
saturates for strong interactions to an almost constant value. 

|| In the open L = 32 system used in the t-DMRG, the density in the middle of the system is n ~ 0.53. 
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We further compare our results to different characteristic velocities of the Bose- 
Hubbard model: (i) the sound velocity in the superfluid regime, (ii) the maximal mode 
velocity of a simplified model in the Mott-insulator, (iii) the maximal mode velocity of 
a fermionic model applicable at low filling. 

(i) For an infinitesimal quench inside the superfluid phase, the rescaled signal 
velocity v s /2 for the spreading of the correlations can be described by the sound velocity 
of the system. However, for finite quenches we expect the rescaled signal velocity to 
be larger than the sound velocity, since modes with higher velocities can be excited. 
In Figure 7 we approximate the sound velocity in the superfluid regime by the sound 
velocity of the corresponding continuous Lieb-Liniger model [36]. For small values of 7, 



it is given by vll = 2n^yy 1 — 2^, where 7 = UJ (2Jn). This expression approximates 
the sound velocity for the Bose-Hubbard model up to an interaction strength 7 < 4 [37] . 

(ii) An idea of the expected signal velocity in the Bose-Hubbard model in the 
Mott-insulating state can be obtained by mapping the system onto a simpler model 
using only three local states, e.g. occupation by n — 1, n , and n + 1 bosons per 
site, where uq is an integer [38, 39]. In the Mott-insulating phase the dispersion 
relation for a particle hole excitations in this effective model can be determined as 
e(Jfe) = ^U 2 - Ue (k)(An + 2) + e (k) 2 /2 [39]. Here e (Jfe) = 2J cos(ka) is the band 
dispersion for the non-interacting case. The group velocity in this case is given by 
v = ■ The maximum velocity v g of a mode given by this model is shown in Figure 7. 
In the strong coupling limit this velocity agrees with the velocity extracted from a 
perturbative calculation where the maximal group velocity is given by aJ/h(2n + 1). 
For the case of n = 1 this results in 3aJ/h. In particular we see that towards the 
critical value in this model the velocity slightly increases. The momentum at which the 
maximum velocity is reached changes compared to the strong coupling limit. 

(iii) For strong interaction and low filling the Bose-Hubbard system can be mapped 
onto a fermionic system [40]. In this fermionic system the velocity is given by 
Vf = 2sin(7rn)(l — 8 Jn cos (nn)/U). 

In Figure 7 we compare our numerical results to the different velocities (i)-(iii). For 
small interaction strength we see that the velocity v s /2 is always larger than the sound 
velocity showing a similar rise with U/J. Comparing further the velocity of different 
quenches in the superfluid regime, e.g. Ui — 1,2 and Uf — 4 in Table 1, the velocity 
v s /2 seems to approach the sound velocity if the parameter changes becomes smaller. 

In the regime of strong interaction the qualitative features of the signal velocity 
are well reproduced by the given approximations. In particular, the velocity is almost 
constant for high interaction values and increases if U/J is lowered towards U/J ~ 6. 
Further the order of magnitude of the values is in good agreement. 

However, the approximation has to be taken with care. It does not take all higher 
energy excitations into account, which will be neccessary to quantitatively describe the 
situation under consideration. Further it does not take into account the initial state, 
i.e. which of the modes are actually excited by the quench. This is in contrast to the 
numerical results in Table 1 which suggest that the value of the observed signal velocity 
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might depend as well on the initial state. In particular if the higher particle fluctuations 
are present in the initial state the observed signal velocity seems to be larger. 

6. Dynamics of entanglement and mutual information 

We turn in the following to the time-evolution of the entanglement and the amount 
of correlations between different subsystems after a quench. As a measure for the 
entanglement we use the von Neumann block entropy (see section 3) and the mutual 
information. Whereas the von Neumann entropy describes the entanglement of a region 
of the system with the remaining part, the mutual information gives a measure about 
the amount of correlations between different subsystems [41]. In the following we first 
analyze the time-evolution of the von Neumann entropy before we turn to the evolution 
of the mutual information. 



U,=2 J, U f =20 J U r 2 J, U f =20 J 

T | , | , | , | i I | , | , | i | i | , | r 




I i I I I i I i I i I i I I i I i I i I i I i I i I 

0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 

timet [J -1 ] timet [J -1 ] 



Figure 8. Time-evolution of the block entropy for different block lengths for a 
quench from Ui = 2J to Uf = 20J. The initial linear rise is the same for different 
block lengths, but depends on the boundary condition of the block (PBC: two interface 
links per block, faster rise of the entropy; OBC: one interface link per block, slower 
rise of the entropy) . For longer times a saturation of the entropy depending linearly on 
the block length is observed. The left panel contains the original time evolution, while 
in the right panel, the obc data has been plotted for times t — 2t, so as to illustrate 
the slower increase of the entropy due to the smaller boundary. The PBC results are 
from ED on L = 14 (circles) and L = 16 (squares) systems, while the OBC data has 
been obtained using t— DMRG up to L = 64. 

von Neumann entropy In Fig. 8 we show the time-evolution of the von Neumann 
entropy for Bose-Hubbard systems with periodic (ED, L = 14, 16 %) and open boundary 

% Due to the computational challenge of calculating density matrices in ED for large blocks we chose 
to limit the local boson occupancy to 3 at most. 
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conditions (t-DMRG, 5 bosons per site are allowed). In both cases we focus on 
contiguous blocks of length I. For open boundary conditions the blocks are aligned 
with one of the boundaries. We plot the difference Sl(1, t) — Sl(1, 0) so that the curves 
for all block sizes start at zero at t = 0. 

In the left panel we show data for ED (/ = 1, 2, 3, 4, 5) and DMRG (/ = 4, 5) results 
for a specific quench from Ui — 2 J to Uf — 20 J. Let us first discuss the ED results: 
Immediately after the quench for t < 0.4J -1 a small dip shows up in all block entropies. 
However, from time t ~ OAh/J up to t*(l) a linear rise of the block entropies can be 
observed. Interestingly all block sizes for I > 1 have the same slope, until they bend over 
to an almost flat behaviour at successively later times. The dip in the entropy at larger 
times is a finite size effect, as can be seen by comparing the data for different system sizes 
(circles for L = 14 and squares for L = 16). The saturation value of the different block 
entropies depends linearly on the block size, defining an entropy propagation velocity 
v e , which is roughly equal to v B determined based on the density-density correlation 
functions in the preceding section 5. In a next step it is instructive to compare the 
entanglement dynamics between different block geometries (two interface links in ED, 
one interface link in t— DMRG). The t— DMRG results are also shown in the left panel 
of Figure 8. The entropy of these open blocks increases more slowly, but converges 
to about the same value at late times as the periodic blocks of the same length. To 
illustrate this convincingly we show in the right panel the same data, but where the 
entropy of the open blocks is shown on a time scale which is twice as fast. Indeed the 
results of the two block geometries agree reasonably well. This nice result lends direct 
support to the picture developed by Calabrese and Cardy [16], where they predicted that 
the saturation of the entropy occurs at t PBC (l) = l/2v for periodic boundary conditions, 
while t* OBC = l/v for blocks aligned with an open boundary. 

These characteristics of the entanglement evolution are similar to the results 
obtained for different models in Refs. [16, 42, 43, 17, 1, 44]. A linear growth of the 
entropy has been seen up to times t = l/2v, where v is the maximal velocity of the 
excitations. Afterwards a saturation of the entropy is seen for t — > oo, and the rate 
of the approach to the saturation value is related to the dispersion relation of the 
underlying model. At a fixed time the entanglement saturates for increasing block 
length, i.e. fullfils a boundary law with the boundary increasing with time. This shows 
that the boundary law for the dynamics of entanglement which has been proven mostly 
for ID spin-systems [17, 18] seems to be valid in more general systems such as the 
Bose-Hubbard system considered here. 

Mutual information While the von Neumann entropy is rapidly increasing in time 
and finally leads to an extensive entropy scaling with the block volume at large times, 
it is interesting to ask whether the vast entanglement entropy also leads to increased 
correlations between subparts of the system. A very useful quantity to address this 
question is the mutual information. The mutual information I(A : B) between two 
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Figure 9. (Color online) Time-evolution of the mutual information (6) between two 
blocks of equal length I — 1 (left panel) and I = 2 (right panel) for different spatial 
separations r. The quench parameters are Ui = 2J, Uf — 20 J, and n = 1. The results 
have been obtained by ED for L = 14 (bold, dashed lines) and L = 16 (thin, straight 
lines) . 



subsystems A and B of the system is defined by 

I(A :B) = S A + S B - S AUB . (6) 

Here Sx is the von Neumann entropy for the subsystem X. The mutual information 
I(A : B) measures the total amount of information of system A about system B [41]. 
Note that AUB is not required to be equal to the total system. Interestingly the mutual 
information can fulfill an area law at finite temperatures, while the entropy is expected 
to follow a volume law [45] + . 

In Figure 9 we show the time evolution of the mutual information between two 
blocks of / = 1 (left panel) or I = 2 (right panel) sites each, shifted by a distance r. 
At t = we expect the mutual information to decay slowly with the block separation 
r, since the starting state is basically a scale invariant critical state and the mutual 
information is just a complicated function of the basic critical correlation functions, 
such as (b\b r ) and (non r ), which all decay algebraically As a function of time t, the 
mutual information is decreasing rapidly with some slow oscillations on top. For later 
times it seems as though the mutual information levels off to a finite value at a time 
which depends again linearly on the distance r. The mutual information at late times 
decays much faster as a function of distance than in the initial state. So even though 
the entanglement entropy is vastly growing in time, this does not lead to enhanced 
entanglement between different subsystems. 

+ In this reference the authors study the case of the total system equal to A U B. 
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7. Conclusion 

In our work we show that correlations and entanglement are very useful quantities to 
characterize the equilibrium and dynamic properties of a quantum many body system. 

In the first part of our work we showed that the static von Neumann entanglement 
signals the quantum phase transition between the superfluid and Mott-insulating state 
without previous knowledge on the type of the phase transition. In the superfluid, the 
von Neuman entropy is in very good agreement with previous predictions by conformal 
field theory [13, 14] with a central charge c = 1. Deviation are only found close to the 
boundaries of the system. These deviations are induced by the inhomogeneous density 
distribution caused by the open boundary conditions. A saturation of the entropy with 
the block length is found in the gapped Mott-insulating phase as predicted [32, 14]. 

In the second part of our work, the time-evolution after a sudden parameter change 
is analyzed with a focus on the spreading of information. Hereby different parameter 
changes are discussed ranging from the change between the superfluid to Mott-insulating 
phase over a quench inside the superfluid to a quench inside the Mott-insulating regime. 
Our study proposes that the Lieb Robinson theorem is valid as well in the considered 
situation of the Bose-Hubbard model. This relies on our findings that a light-cone 
like evolution takes place in different correlation functions. The velocity of the front 
evolving in the correlation functions is discussed and compared to different characteristic 
velocities of the system. The validity of the Lieb-Robinson theorem is further supported 
by the von Neumann entropy which shows for a certain block length a linear growth 
for short times. To a good approximation the rate of the block entropy increase is 
shown to follow a boundary law, while the entropy value at saturation seems only to 
depend on the block volume. Our results nicely corroborate recent analytical results 
[16, 17, 18]. However, in contrast to the von Neumann entropy after a quench from the 
superfluid to the Mott-insulating regime, the mutual information between two spatially 
separated small blocks relaxes to a lower value than in the starting state. So even 
though the entanglement entropy is vastly growing in time, this does not necessarily 
lead to enhanced entanglement between different regions of the system. 
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